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We calculate the spectrum of quasiparticle excitations in the core of isolated pancake vortices 
in clean layered superconductors. We show that both the circular current around the vortex center 
as well as any transport current through the vortex core is carried by localized states bound to 
' the core by Andreev scattering. Hence the physical properties of the core are governed in clean 

| high-K superconductors (e.g. the cuprate superconductors) by the Andreev bound states, and not 

y—i . by normal electrons as it is the case for traditional (dirty) high-K superconductors. 

a ' 
3 ■ 

■ I. INTRODUCTION 

m : 

' We discuss specific aspects of the core of a vortex line in layered high T c superconductors. The physics of these 
vortices is governed by two distinct length scales, the London penetration depth in the planes, An ~ 10 3 A, and 
the coherence length in the planes, £|| 10 — 20A. The penetration depth is the electromagnetic length scale of a 
. vortex. The physics on this length scale is well described by a combination of macroscopic electromagnetism, London's 
\Q ' theory for supercurrents along the layers, and interlayer Josephson coupling. This description breaks down in the 
r—i , core of the vortex, i.e. at distances of order £y from the center of the vortex. Thus, physical properties of the core 
' carry information on the microscopic physics of high T c superconductivity. The small coherence length of high T c 
superconductors makes the vortex core a good potential sensor for microscopic mechanisms of superconductivity. Our 
discussion of the vortex core in high T c superconductors is based on the Fermi- liquid model of superconductivity. The 
physical properties of the vortex core predicted by this model are spectacular, unique, and could serve as fingerprints 



On 



, of the traditional pairing theory of superconductivity. _ 

S The vortex core of traditional high-K superconductors is well described by the Bardeen-Stephen modelEJ which 
represents the core by a region of normal electrons. The Bardeen-Stephen model is justified as long as the mean free 
■ path, £, is much shorter than the core size, so that the motion of an electron gets randomized before it leaves the 
core. This condition is not fulfilled in high T c superconductors which are generally clean superconductors with I > £y . 
The coMjjjpf a vortex in a clean superconductor was first studied in the classic papers of Caroli, Matricon, and de 
GennesHtil These authors calculated the spectrum of quasiparticle states in the core, and showed that electrons and 
^ . holes form bound states at energies below the bulk energy gap. Further early studies of the excitation spectrum in 
the core can be found in Refs. 

More recent theoretical work was stimulated by the direct,p.bservation of core states in NbSe2 by scanning tunnelling 
spectroscopy (STS).Bi3 The recent report of STS in YBCOcM provides new information on the excitation spectrum of 
vortices in the high T c cuprates. Consequently-theoretical efforts focused on the tunnelling density of states of bound 
states in isolated vortices and vortex lattices. E9TL3 These calculations show that the bound states in the core have a 
different nature compared with the usual quantum mechanical bound states in a potential well. The core states are 
coherent superpositions of particle states and hole states and are formed by repeated Andreev scattering from the 
pair potential (order parameter) in the core. Andreev scattering is a process of "retroreflection" of excitations: spatial 
variations of the amplitude or the phase of the order parameter induce branch conversion of electron-like excitations 
into hole-like excitations, and vice versa. Bound states occur at energies at which the phases of multiply reflected 
electron-like and hole-like states interfere constructively. The charge current carried by an incoming electron and an 
outgoing Andreev reflected hole is identical because the reversal of the velocity in an Andreev reflection process is 
compensated by the reversal of the charge due to electron-hole conversion. Consequently, Andreev bound states can 
transport a charge current, unlike bound states in a potential well. Charge conservation requires that the current 
carried by the bound states inside the core is transported outside the core by bulk supercurrents. This leads to an 
interplay between supercurrents flowing past the core and the bound states in the core. Hence, the physics of the 
'normal core' in clean superconductors is basically the physics of the bound states in contact and intimate exchange 
with the superconducting environment outside the core. 

Consider a stack of "pancake" vortices forming an isolated vortex line whose axis is oriented perpendicular to 
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the layers. We investigate the current distribution in the core of a pancake vortex, and show how this distribution 
changes if the vortex is exposed to a bulk supercurrent, or the circulation is changed from 2ir to An. We calculate the 
spectral current density, which carries the information on the contribution of the states in a given energy interval to 
the total current density. A supercurrent in homogeneous superconductors is distributed over all continuum states. 
These states exhibit Doppler shifts of their energies, Se — Vf ■ p s , in the presence of a phase gradient in the order 
parameter (or superffuid momentum, p s = §V%— - A). The total current is obtained by adding the contributions of 
states with positive shifts from quasiparticles co-moving with the flow and the contributions with negative shifts from 
quasiparticles that are counter-moving relative to the flow field. We find that the currents in the core have a very 
different spectral distribution from bulk supercurrents. The continuum states (scattering states) show smeared out 
Doppler shifts, and contribute very little to the total current. The dominant contributions to the circulating currents 
around the vortex center, as well as the currents through the core, come from Andreev bound states. Hence, the 
physics of vortex cores in clean superconductors (£y <C £) is very different from the physics of the vortex core in a 
dirty superconductor (I <C £11), which is well described by a continuum of normal electronic states. The calculations 
presented in this paper concentrate on stationary properties of the vortex core of clean layered superconductors. We 
expect more spectacular effects in the dynamic properties. The bound states react sensitively to the environment 
outside of the core. This leads to a coupling of the collective degrees of freedom in the London range of the vortex 
and the bound states in the core, which will produce a rich spectrum of largely unexplored dynamical phenomena. 

Below we present analytical and numerical calculations for the states in the vortex core. We use two versions 
of a quasiclassical formulation of the BCS theory— pi superconductivity: a) Andreev's theorjO wiiich represents 
the quasiclassical limit of Bogolyubov's equations,ll3 and b) the quasiclassical theory of EilenbergeiO, Larkin, and 
Ovchinnikovo which represents the quasiclassical limit of Gorkov's Green's function theory. Andreev's theory and 
the quasiclassical theory are essentially equivalent for clean superconductors, and in this limit the choice of approach 
is largely a matter of taste. However, the quasiclassical theory has a wider range of application. It is the general- 
ization of Landau's Fermi liquid theory to the superconducting state, and is capable of describing a broader range 
of superconducting materials and phenomena, such as dirty superconductors or superconductors with short inelas- 
tic lifetimes (strong-coupling superconductors). Ell Section II contains analytical results for the bound states and the 
spectral current density_£or a pancake vortex with a superimposed bulk supercurrent. These results are obtained from 
Andreev's HamiltonianO by the methods described in Ref. |22|. In section III we discuss the numerical results, which 
are obtained using the quasiclassical theory of Fermi-liquid superconductivity. We solve the quasiclassical transport 
equations to obtain self-consistently the pair amplitude (order parameter) for pancake vortices. Given the pair ampli- 
tude we calculate the excitation spectrum in the core of the vortex, and deduce from it the spectral current density. 
The numerical calculations are done for layered superconductors with s-wave pairing. The analytical and numerical 
results confirm and complement each other, and they establish the important role of the bound states for the currents 
in the core region of a vortex. 



II. ANALYTICAL RESULTS 



In this section, we investigate the spectrum of current carrying states of a two-dimensional pancake vortex in 
equilibrium at temperature T, in the_presence of an externally imposed supercurrent. We ignore the spin degree of 
freedom of a quasiparticle excitation;Ej in this case it is sufficient to work in the two-dimensional space of particle-hole 
degrees of freedom. Operators in this space are 2x2 matrices, and we use the notation fi, f%, T3, for the three Pauli 
matrices hvnarticle-hole space (Nambu space), and 1 for the unit matrix. The Hamiltonian for the quasiparticle 
excitations,E-3 the Bogolyubov Hamiltonian, then reads: 

Hb = MP + -A(r))f 3 + A(r), h (p) = ( ^— ^ J , (1) 
c \ 2m / 

where p = (p x ,p y ) and f = (x,y) are the momentum and position operators appropriate to a particle moving in two 
dimensions, pf = mvf is the Fermi momentum and A is the electromagnetic vector potential. The order parameter, 
A(r) is an off-diagonal matrix and is generally represented by a linear combination of fi and t-i- 

In the absence of an externally imposed supercurrent, we write the order parameter of the vortex as 

A (r) = A F(r)fi exp(ip f 3 ) , (2) 

where Ao is the magnitude of the order-parameter of a bulk superconductor at temperature T, F(r) is the normalized 
profile of the vortex, which is a monotonically increasing function of r obeying F(0) = 0, F(oo) = 1, and ip is the 
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angular coordinate of r with x = r cos tp, y = r sin ip. The main assumption we make in this section is that the 
order-parameter in the presence of a superfiow p s = h/2 V% — e/c A has the form 

a i * % 

A(r) = exp(+-Vx ■ rf 3 )A (r) exp(--Vx ■ rf 3 ) (3) 

We assume throughout this section that p s is small compared to the bulk critical current, Vfp s <C Ao. 

The principal physical quantities with which we shall concern ourselves are the spectral current density, and the 
total equilibrium current density which is related to j(r, e) by 

j(r,T)= [ de](v, e)f(e), f(e) = 1 ■ (4) 

J exp(e/T) + 1 

We shall also make reference to the local density of states, N(r, e). The quantities j(r, e) and N(r, e) may be expressed 
in terms of the one-particle Greens function, l/(e — Hb) or, equivalently, the "spectral function" 5(e — Hb). Using 
the spectral function, we find that in Dirac notation, 

j(r, e) = 2e<r|{P±|^, S(e - H B )}\r) 1A , (5) 

N(v,e)=2(r\S(e-H B )\r) hl , (6) 

where the subscript 1, 1 denotes the upper left element of the 2x2 matrices, thereby selecting out the particle sector 
of the spectral function and the factor 2 takes into account both spin projections of the quasiparticles. 



A. Andreev Hamiltonian 



Most calculations of the properties of superconductors with inhomogeneous order parameters are simpler in the 
quasiclassical limit, where one takes advantage of the separation in the scales of the wavelength of quasiparticles 
near the Fermi energy and the characteristic scale for spatial variations of the pair potential, i.e. h/pf <C £o- The 
quasiclassical limit of the Bogolyubov Hamiltonian ([!]) is the Andreev Hamiltonian in which the kinetic energy in 
(lj) is replaced by an operator that is linear in the gradient £ZI Let us define the normal-state density of states at the 
Fermi level, Nf = pf/2irvf, and introduce the directions, k = (cos^t,, sin and 1 = (— sin tpi-, cos <Pk), which are, 
respectively, parallel and perpendicular to trajectories of a quasiparticle wavepacket in the quasiclassical description, 
i.e. Vf — f/k. The coordinates along these directions are defined by r = £k + 77I. In addition we work in the limit 
A/£ 3> 1, in which case the vector potential is approximately constant in the vicinity of the vortex core and can be 
neglected. The order parameter in (|^) can be written as 



A(r) 



exp 



-ip s • kCf 3 



(7) 



By performing a gauge transformation that removes the factor of U in Eq. (|7j) we obtain the spectral current density 
and the local density of states in terms of the Andreev Hamiltonian for an isolated vortex, 



F(f) 

v fP(T 3 + A — — (fiC + T 2 T]), 



(8) 



j(r, e) ~ Airev'jNf 



2tt 



Zir 



[H A + v fPs -k})\0 



1,1 



(9) 



N(r,e) ~ Airv/Nf 



Z7T 



[H A + v fPs -k})\Oi,i 



(10) 



where |C) is an eigenvector of the "one-dimensional" trajectory coordinate operator, (: Q() — CIC)- The operators C 
and k • p = p^ appearing in Ha are canonically conjugate: [CiPcl = The quasiclassical interpretation given to (|^) 
is as follows: quantum-mechanical evolution in particle-hole space takes place along classical trajectories parallel to 
k having a fixed value of 77 = 1 • r. Thus, rj is identified as a c-number impact parameter £3 
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B. The current density of a vortex in a flow field 



Let us write the current density at temperature T as 

/oo 
rfej(e,r)/(e), (11) 
-A 

where A is a high energy cutoff that serves to make manipulations of j(r.T) well defined; large positive energies 
are automatically cut off by the Fermi function, /(e). Where no ambiguity arises, we shall take A = oo. Defining 
q = w/Ps, and using (||) we have 

j(r.T)=4irev 2 f N f ^± T dek((\S(e - H A )\()i,if(e + q-k) . (12) 

Next, we split up the energy integrals into the following terms 

j(r,T) = j!(r,T) + j 2 (r,T) + j 3 (r,T) , (13) 

jx(r,T) = AnevjNf ^ ^k (q-k) (C\S(-A H A )\Qi,i , (14) 
j 2 (r,T) = 4nevjN f ^ J°° dek{C,\5{e - H A )\Qi,i [/(e + q-k) - /(e)] , (15) 

j 3 (r,T) = 4nev}N f J J dek{C\5(e - H A )\() h if(e) . (16) 

The three contributions to the current have different interpretations. 

1. Since A is large, (£|(5(— A — H A )\C)i,i may be replaced by its high energy, normal-state limit, \j1-KVj and 

ji(r,T) = e^JV>p s . (17) 

This term coincides with the T = current of a uniform superconductor. 

2. The term-^r,! 1 ) contributes to "backflow", since it always yields a current with a component in the — p s 
direction.L3 The term j2(r,T) contains the current carried by the bound states and also a correction to the 
T = current due to the thermal breaking of pairs. To appreciate these points we look at this term in two 
limits, assuming q <C Ao. (i) For T = we have 

j 2 (r,0)=4 7 re^iV / F ^ [ dek {C\S(e - H A )\Q hl . (18) 
Jo zn Jo 

The small size of q ensures that the energy integral only selects states in the gap, thus j2(r,0) only obtains 
contributions from the bound states, (ii) For T / and assuming Ao(r) to be that of a uniform system, 
Ao(r) = A fi, we write /(e + q-k) — /(e) w q-k/'(e) and obtain 

j 2 (r,T) = -evjNjPsJ^de^M^e^ - A 2 )[-f'(e)} = -evjN fPs Y((3A ) , (19) 

where Y(f3 Aq) is the Yosida function which gives a quantitative measure of the thermal breaking of Cooper 
pairs. 

3. The term j 3 (r,T) is independent of p s and is simply the current density of a vortex in the absence of an externally 
imposed supercurrent. 
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C. Current carrying bound states at the center of the vortex 



Consider the spectral properties at the center of the vortex. For rj = 



r;=0 



vfPch + A F(()ti 



(20) 



with F(— C) = — F(C) accounting for the tt phase change across the vortex core. This special case of the Andreev 
Hamiltonian is identical in form to the continuum Hamiltonian used to describe irans-polyacetylene containing a single 
topological soliton.EH It is known that this Hamiltonian always has a non-degenerate bound state at zero energy.ca 
Whether or not it has other bound states depends on the form of the profile, F((). For the single quantum vortex 
and trajectories through the center there are no other bound states. The eigenfunction for the zero-energy bound 



state, "00 (C) i is found by solving 



AoF(Qti ipo(C) = 0. The normalized solution is 



MO 



exp 



Vf Jo 



d('F(C 




(21) 



e?C exp 



2A 



v f Jo 



dC'F(C) 



(22) 



where L is a profile dependent quantity with the dimensions of leng th, L ~ vj / Aq. Analytical estimates of the bound 
states at distances far from the vortex are given in appendix IV B. 

For energies |e| < A only the bound state of Ha will contribute to the spectral current density (and the local 
density of states) , 
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dip k 



j(e,0) ~ inevjNf \ k<J(e - v fPs ■ k) [^,(0)^3(0) 



i,i 



(23) 



e 2eN f A Q((v fPs ) 2 -e 2 ) 
Ps L/£ ^(v f p s ) 2 - e 2 



P s , |e| < Ao • 



(24) 



There is a simple relation between j(e, 0) and N(e, 0) when |e| < Ao. In (p3]), the delta function in the integrand of 
j(e, 0) effectively replaces k by (e/vfp s )p s . Taking this factor outside the integral leaves an integral identical to that 
of the local density of states. Consequently, 



j(e,0) = e— iV(e,0)p a , |e| < A . 

Ps 



(25) 



Note that the contribution of negative energy (bound) states to the total current density lies in the — p s direction, 
i.e. opposite to the externally imposed supercurrent. At zero temperature the total current density originates from 
the bound states having energies — A < e < 0, 



2eiV f A 

jbound(0,T = 0) = / j(e,0)de= 777^v f p s 

'-Ao L IK 



(26) 



The current density of an isolated vortex with p s — vanishes at the center of the vortex, i.e. j3(0,T) = 0. We can 
combine the result in (|2^) with ji given in ( [Tt| ) to obtain the total current density at T = 0: 



j(r = 0,T = 0) = e^A^p 



2A /e A 



-VfPs 



(27) 



Thus, for sufficiently small p s the bound-state contribution dominates ( |27| ) and j(r = 0,T = 0) will point in the — p s 
direction. 
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D. Particle conservation 



For equilibrium conditions the divergence of the current density vanishes. From (17), ji has a vanishing divergence, 



and for an undisturbed vortex we have V • j3(r) = 0. Thus, S(r) = V • j(r) = V • J2(r). At T = 0, 



S(T)=4nev f Nf I q dev f -^({(\6(e-H A )\0i 



In Eq. (49) of appendix IV A we show that 



d - F(r) 

v f —{<;\8(e-H A )\Q lll = A ^ r ±i X 



yielding 



S(r) = AnevjAoNf 



2tt 



(Wi-irMKI^-ftOlO 



Since q <C Ap, only the bound state, 4>o(('iV)i contributes to the e integral. Thus, 



2tt 



S(t) = iwev f N f Ao / ^/ deS(e - e (r?)) r,)^-{U ~ hv)MC, v) : 



qk 



F(r) 



(28) 



(29) 



(30) 



(31) 



where eo(??) is the bound-state energy for an impact parameter 77. At small distances from the center of the vortex 
(r <C £) F(r) « r/^, with ^ w Vf/2A , and the lowest energy bound state is then 



1 1 



exp 



(32) 



eo(»?) 



A V 
-A 0? . 



(33) 



Substituting -0q and £q into ( pl|) and setting exp(— ^C 2 /^ 2 ) sa 1 yields 



S(r) 



q 

A 



(34) 



which is non-zero, indicating that the ansatz (|^) is not physically correct for any value of p s . This failure to satisfy 
the conservation law is due to the lack of self-consistency of the order parameter in Eq. (||) in the presence of the flow 
field. In the absence of pinning, the vortex will move in response to a flow field, even one of arbitrarily small strength. 
The results in section II implicitly assume a pinned vortex. Thus, there will be distortion of the vortex away from its 
cylindrically symmetric equilibrium form (||). In the following section we show that the self-consistently determined 
vortex order parameter, which includes the deformation by the flow field, restores the conservation law. 



III. QUASICLASSICAL RESULTS 

A versatile and efficient nictiipd|jc)r calculating local spectral properties of superconductors is the quasiclassical 
theory of superconductivity.t3E3cjH This theory is well adapted for a numerical approach to microscopic problems 
in superconductivity, such as the calculation of the structure and the excitation spectrum of vortex cores for su- 
perconductors with isotropic, anisotropic or unconventional order parameters. The quasiclassical theory is the only 
theoretical formulation which can handle equally well clean and dirty superconductors, as well as more complicated 
geometries than that of an isolated vortex with cylindrical symmetry or a perfect vortex lattice. It can be interpreted 
as the generalization of Landau's theory of normal Fermi liquids to the superconducting state. The quasiclassical 
theory shares with Landau's theory the semiclassical description of the orbital degrees of freedom of quasiparticle 
excitations. On the other hand, the internal degrees of freedom, i.e. the spin and the particle- hole degrees of freedom, 
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are described by quantum mechanics. Quantum mechanical coherence of particle and hole excitations is the basis of 
all superconducting effects such as persistent supercurrents, flux quantization, Josephson effects and Andreev reflec- 
tion. Here we use the quasiclassical theory for our investigations of the vortex core. Numerical work on vortices in 
superconductors using the quasiclassical theory started with a series of publications by Kramer, Pesch, and Watts- 
TobinflH'tHl More recent work includcfi-pinmng of vortices at small defects,c3 vortices in .superfluid 3 He and other 
systems with unconventional pairing ftltB the excitation spectrum of bound quasiparticles,E3E3 and the spectrum of 
moving pancake vortices E3 

We use in this article the notation of Refs. pqj39| , pi]. The central objects of the quasiclassical theory of supercon- 
ductors in equilibrium are the quasiclassical propagators g R ' A {pf, r; e), which are 2x2 matrices in the particle- hole 
index, 

±R,A ( 9 R ' A (Pf,r;e) f R > A (p f ,r;e)\ 



r--( P/ ,r; e ) £^( P/ ,r;e) 

The variables e and p f stand for the energy of an excitation and its momentum (on the Fermi surface). The momentum 
variable reduces to pf = pf(cos(pk, sin^fc) for an isotropic Fermi surface in two dimensions (see section If). General 
symmetries lead to the following fundamental relation between g R and g A , 

<1 X -;(.,/•') T r ;i . (36) 

We use, as described in Section II, the notation fi, t<i, T3, for the three Pauli matrices in particle-hole space, and 1 
for the unit matrix. The off-diagonal terms f RA in Eq. ( |35| ) arc the pair amplitudes. They vanish in the normal 
state, and measure the amount of particle-hole mixing in the superconducting state. The diagonal elements of the 
propagators determine the density of states, 

N( Pf , r; e) = Nf 9*<Pf,r^-£(p f ,rrf (3?) 

and the equilibrium current density. The most detailed information on the current distribution is obtained from the 
spectral current density, 

j{pf,r;e) = ev f Nf (A/+(p/,r; e) - W- (p/, r; e)) , (38) 

where A/±(p/,r;e) = A(±p/, r; e)/iV/ is the dimensionless density of states for co-moving (+) and counter-moving 
(— ) excitations along the trajectory line defined by p/, and v/ is the Fermi velocity at the point pj on the Fermi 
surface. This spectral density measures the contributions of quasiparticle states with energy e and momentum near 
the Fermi surface point p/ to the current density at position r. The full current density is obtained by weighting 
the spectrally resolved current density by the occupation probability of the quasiparticle states, then integrating over 
Fermi momenta and energies. For equilibrium states, 

j(r) =2jdejdp f j( P/ , r; e) (/(e) - , (39) 

where /(e) is the Fermi distribution function. The symbol J dpf denotes a weighted integral over the Fermi surface. 
The weight at p/ is cx| v/ | , and the integral is normalized, / dpf 1 = 1. The spectral current density is particularly 
well suited for our study of the importance of Andreev bound states for the current flow in a vortex core. These 
bound states appear as delta functions in the spectral current density at energies below the bulk energy gap. The 
spectral weight of the delta function, combined with the occupation of the bound state, determine its contribution to 
the total current density. _ 
We calculate <7 fl ( P /,r; e) from Eilenberger's transport equationlij 

(e + -w { • A(r))f 3 - A( P/ , r), g R (p f , r; e)] + thv f ■ Vg R {p fl r: e) = , (40) 

supplemented by the condition of analyticity in the upper half of the complex e-plane, and the normalization condition 

g«( P/ ,r;e) 2 = -7r 2 i. (41) 

For a fixed Fermi momentum pf this is a first order differential equation along a straight-line classical trajectory in 
the direction of the Fermi velocity v/. The propagator g R (pf, r; e) at a chosen point of interest, r, is determined by 
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the solution of ( |4C| ) along the trajectory through r in the direction v/. Complete information on the local physical 
properties at point r, such as the current density, is obtained by sampling all trajectories through r. The propagator 
g R (pf, r; e) is intimately related to the 2x2 density matrix of the particle- hole degrees of freedom of a quasiparticle 
moving along the classical trajectory specified by p/, r. Thus, g R (pf, r; e) describes the state of the internal degrees of 
freedom of the excitation. The internal state, i.e. the amount of particle- hole mixing, may change along the trajectory 
as a consequence of the off-diagonal pair potential, A(p/,r), which acts as a driving term that 'rotates' the particle- 
hole pseudo spin. The pair potential couples particle and hole excitations, and is the origin of particle-hole coherence. 
It depends on the real space position, r, and, for anisotropic superconductors, on the Fermi surface position, p/, 

^>=(-A.(P, r ) A(P 0' r) )' (42 » 
The pair potential must be calculated self-consistently from the gap equation, 

A( P/ ,r) = J dp/y(p/,p/) J ^lm/ fl (p/,r;e)(l-2/( e )), (43) 

where V(p/, p^) is the pairing interaction, which determines the orbital symmetry of the pair potential, its magnitude 
and T c . 

Our procedure for numerical calculation of the currents in the core of 2D pancake vortices is the following. We 
first solve self-consistently the gap equation and Eilenberger's equation at Matsubara energies. This allows us to 
determine the pair potential and the supercurrent density. We then insert the pair potential into Eilenberger's 
differential equation at real energies, and obtain from its solution the excitation spectrum: the density of states and 
the spectral current density. The differential equations are solved by standard 4 th order Runge-Kutta and multiple 
shooting methods, and self-iconsistency is achieved iteratively by using alternatively a relaxation method and the 
Mobius-Eschrig algorithm.^! 

We consider three examples of pancake vortices: isolated, i) singly-quantized and ii) doubly-quantized s-wave 
vortices, and iii) a pinned s-wave vortex in the presence of a uniform transport supercurrent. We choose a temperature 
of T — 0.4 T c , unless otherwise noted, and assume k = A/£ ^> 1. In this limit the vector potential is essentially constant 
in the core region, and can be neglected. 



A. Spectral current density of a singly-quantized s-wave vortex 

Figure 1 shows the amplitude of the order parameter of a singly-quantized s-wave vortex. The amplitude is isotropic 
and vanishes linearly in the core. The variation of the amplitude and phase along two trajectories are also shown in 
Fig. 1. For trajectory a passing through the center of the vortex, the phase changes discontinuously and the amplitude 
vanishes linearly at the vortex center. For trajectory b, with impact parameter 77 = 3.0£o, there is only a small change 
in the amplitude of A. For singly-quantized vortices the phase of the order parameter is the more important factor 
determining the spectrum of bound states. 

Figure 2 shows the angle-resolved local density of states for the two trajectories shown in Fig. 1. For trajectory (a) 
through the center of the vortex, the spectrum shows a zero-energy bound state separated from the continuum that 
begins at the bulk gap. The bound state results from constructive interference of particle- and hole-like quasiparticles 
that undergo Andreev reflections from the vortex order parameter .-.This bound state corresponds to the zero angular 
momentum bound state found by Caroli, de Gennes and MatriconH'fj A zero-energy bound state is always present for 
trajectories in which the order parameter is real (up to a constant phase factor) and has different signs when going 
to ±00 along the trajectoryEa 

Bound states with non-zero energies are found for trajectories with a nonzero impact parameter measured from 
the vortex center. Theserbound states correspond to the spectrum of bound states with non-zero angular momenta 
obtained by Caroli et alfl Figure 2b shows the spectrum for a trajectory with an impact parameter of r\ = 4.2 £0 
and ■ p s (r) > measured at the point of closest approach to the vortex center. The bound state is shifted down 
in energy to e/2irT c ~ —0.22, and the continuum states are shifted and inhomogeneously broadened by the Doppler 
energy, Ae = vt ■ p s (r). The spectrum near the onset at point 1 in Fig. 2b has low weight and corresponds to the 
continuum edge at e = A far from the impact point, while the peak in the spectrum at point 2 corresponds to the 
maximum Doppler shift, e = A + v/ -p s (R) at the impact point R. Note the development of the BCS coherence peak 
as the density of states is sampled further from the vortex center. 
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FIG. 1. The magnitude of the pair potential, [A(r)|/2-7rT c , at T — 0AT c for a singly quantized vortex 
in an s-wave superconductor is shown in the 3D plot. The 2D plots show the order parameter amplitude 
and phase of the order parameter along a trajectory (a) passing through the center of the vortex, and (b) 
along a trajectory with an impact parameter of rj = 3.0£o- The order parameter is real along trajectory (a) 
and the phase changes discontinuously by n. Along trajectory (b) there is little change in amplitude, but a 
substantial, continuous change of phase. 



The density of states of an s-wave vortices has been investigated by several authors Our emphasis is on 
the importance of the Andreev bound states for the current distribution in the vortex core. We show in Fig. 2c the 
spectral current density for the trajectory with 77 = 4.2 and Vf ■ p s (r) > 0. The net current carried by the states 
at the point ±p/ on the Fermi surface is obtained by weighting this spectrum by the equilibrium distribution and 
integrating over all energies. Thus, for T — > only the negative energy states contribute. It is clear from Fig. 2c that 
the current in the vicinity of the vortex core is carried almost entirely by the bound states with — |A| < e < 0. The 
continuum states give almost no net contribution to the current in the core. Figure 2d shows the spectral current 
density of the set of bound states with trajectories v/ = ±«/y as a function of the impact parameter 77 for < 77 < 6£o . 
The peak at e/2nT c ~ —0.027 corresponds to the trajectory with impact parameter 77 = 0.2£o- The bound state energy 
decreases with increasing distance from the core. For small 77 we obtain, €0(77) ~ — 2(r]/£o)A, in reasonable agreement 







with the analytic estimate in Eq. (J33|) . As indicated in Fig. 2d the contribution of the bound state to the current 
density decreases as the impact parameter increases. However, even at a relatively large distance, r\ = 6£o, the bound 
state still contributes significantly to the circulating current density of the vortex. 

The evolution of the bound state energy for small impact parameters can be written in terms of the angular 
momentum of an excitation about the vortex center, C. z — Pfr/; i.e. ea(jl) = —C z u>o, where %uiq = 2?iA/p/£o <C A. 
This spectrum was originally obtained by Caroli et al.H by solving the Bogolyubov equations. In the Bogolyubov or 
Gor'kov formulation the spectrum is discrete: C z = (m + |) K with m = integer and Hojo defining the level spacing 
of the low-lying bound states in the core. The lowest energy bound state in the core has a zero-point energy of 
e = ~ A 2 /_E/ <C A which is outside the resolution of the quasiclassical or the Andreev theory. The discrete 

spectrum of the Bogolyubov theory corresponds to the continuous Andreev spectrum in the limit where the level 
spacing is small compared to all other relevant energy scales, i.e. Hujq <C fcsT , h/r, etc. This is generally an 
excellent approximation in conventional type II superconductors. For the high T c cuprates the discrete level structure 
is expected to play a more important role, particularly— in the transverse response of vortices in the ultra-clean limit, 
i.e. ljq ^ 1/t, where r is the mean scattering timecicHl 
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FIG. 2. a) Local density of states at the center of a vortex for a trajectory passing through the center of 
the core. The width of the bound state is set at 7/2-71- T c = 0.0004, the continuum edge is at e = ±A, and 
the temperature is T — 0AT c . b) Local density of states at R = (4.2, 0)£o for the trajectory V/ = (0, 
The bound state is shifted, e/2iiT c ~ —0.22, and the continuum states show the Doppler broadening, c) The 
spectral current density for the same position and direction as in b. The Fermi function for T = 0.47^ is 
also shown. Note that the current density is dominated by the negative energy bound state, d) The spectral 
current density for a set of parallel trajectories as a function of impact parameter for < r\ < 6£o- The 
spatial separation between neighboring trajectories is 0.2£o- 
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B. Spectrum of a doubly-quantized s-wave vortex 



It is interesting to compare the single-quantum vortex with the axially symmetric, 4n vortex, A(r) = | A(r)| exp 2i<p. 
The double quantum vortex has higher energy than a pair of isolated single-quantum vortices; however, once created 
the double-quantum vortex is metastable against dissociation into singly-quantized vortices. The amplitude of the 
order parameter for the double-quantum vortex decreases as |A(r)| <~ r 2 for r < £o as shown in Fig. 3a. 




FIG. 3. a) The amplitude of the order parameter for a 4tv vortex at T = 0.4T C . Note the quadratic 
behavior for r <C £o- b) Local density of states at the center of the vortex for a trajectory passing through 
the center of the core. Two bound states are present at energies, e/2irT c = ±0.18. c) The plot of J y (x,0) vs. 
x shows a reversal of the direction of the current for x < 1.9£o. d) The magnitude of the current density for 
the 47r vortex. The corresponding quantities for the 2-k vortex are shown for comparison (dotted curves). 



In contrast to the 2-7T vortex there is no sign change of the order parameter for trajectories passing through the 
center of the vortex core. This difference has a profound effect on the spectrum of Andreev bound states in the core. 
Fig. 3b shows the excitation spectrum of the doubly-quantized vortex at the center of a trajectory passing through 
the center of the vortex core. A symmetric spectrum of two bound states at e±/2irT c = ±0.18 are separated from the 
continuum. Figure 3 also shows the current density of the doubly-quantized vortex. The remarkable feature is the 
reversal of the current direction in the core, i.e. for r < 2£o (see Fig. 3c and 3d). This current anomaly is associated 
with the appearance of a counter-moving Andreev bound state below the Fermi level (e = 0). The evolution of the 
spectral current density is shown in Fig. 4. The trajectories are parallel to y and the spectral current density is 
shown as a function of the impact parameter. At distances greater than x ~ 2£o two bound states lie below zero 
energy, and both states are co- moving with the circulating phase gradient, p s . As the vortex core is approached the 
co-moving bound state nearest the Fermi level moves to higher energy, and a counter-moving bound state above the 
Fermi level (not shown) moves to lower energy. These two states cross the Fermi energy (e = 0) at approximately 
x = 2£o, leading to a reversal of the integrated current density inside the core. The cumulative current density for 
each trajectory is shown as the thick solid line in each panel of Fig. 4. 
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FIG. 4. The spectral current density of the 4n vortex for impact parameters, x = 0.5£o, ...,3.5£o- The 
cumulative spectral weight is shown as the thick solid line in each panel. Note the appearance of the 
counter-moving bound state at x = 2.0£o and the corresponding reversal in the integrated spectral weight 
for x < 2£o- 



C. Spectrum of a pinned s-wave vortex in a transport current 

Finally, consider the current and excitation spectrum of a 2ir vortex in the presence of a uniform supercurrent 
jtr = jtr-x. In the absence of pinning the vortex will move in the direction (— y) in order to reduce the kinetic energy. 
Thus, in order to investigate the excitation spectrum in the presence of a transport current we must pin the vortex to 
the lattice. Our model for the pinning center is a normal metal inclusion where the pairing interaction (or the local 
T c ) vanishes. 
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FIG. 5. a) The amplitude of the order parameter for a pinned 2ir vortex at T = 0.4T C . The normal 
inclusion has a diameter of 0.4£o, and the imposed transport current corresponds to p s = 0.025 ~ 1 (1, 0). The 
dotted curve corresponds to | A(ar) | in the absence of the normal inclusion, b) Current density of the pinned 
vortex for a trajectory passing through the center of the core, c) The bound state spectrum at the center of 
the core of the pinned vortex in a uniform flow field. The two nearly zero-energy bound states correspond 
to co-moving and counter-moving trajectories, d) The spectral current density for a trajectory through the 
core. The negative energy counter-moving bound state carries the backflow current in the core. The thick 
line is the cumulative spectral weight for the current density. Note the scale changes for the two panels. 



Fig. 5a shows the order parameter of a pinned vortex for an s-wave superconductor and a pinning center with a 
radius of 0.4£ - In the presence of a transport current the amplitude of the order parameter deforms; it is suppressed on 
the high current side of the vortex as shown in Fig. 5a. For relatively small transport currents, e.g. p s = 0.02^ 1 (1, 0), 
the center for the phase winding lies within the normal inclusion. However, as shown in Fig. 5b, the vortex current no 
longer vanishes at the center of the vortex core; there is substantial current through the vortex core region, including 
the normal inclusion. The current density inside the normal inclusion is carried by the Andreev bound states and is 
a consequence of the proximity effect. The bound state spectrum at the center of the vortex is shown for a trajectory 
parallel to the transport current and passing through the vortex center. The negative energy bound state carries the 
transport current inside the normal inclusion. Fig. 5d shows the spectral current density measured at the center of 
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the normal inclusion for the trajectories with v/ || p s . Note that the bound state dominates the current and that 
this current is opposite to the applied transport current. This result was also obtained in section H, without taking 
into account the distortion of the vortex core order parameter. This led to a violation of charge conservation in the 
core. Our numerical calculation shows that the main features of the analytic model for the bound state spectrum and 
the self-consistent determination of the order parameter for the pinned vortex in the presence of a transport current 
guarantees that charge is conserved. 

CONCLUSIONS 

We have discussed the current carried by the excitations of s-wave vortices in clean layered superconductors. The 
spectral current density was introduced in order to identify the excitations that determine the transport and circulating 
currents of a vortex. The bound states of the vortex carry most of the current in the vicinity of the core, including 
transport currents that flow through the core of a pinned vortex. Far from the vortex core currents are carried 
primarily by the bound-pair continuum that forms the condensate. For currents flowing through a pinned vortex, 
current conservation is maintained by "spectral transfer" of the current carried by the Andreev bound states to the 
continuum states outside the core. A novel example of the evolution of the spectral current density is provided by 
the double quantum vortex which shows the connection between the spectrum of bound states and the symmetry or 
topology of the order parameter. At low temperatures (T = 0.4T C ) the double quantum vortex exhibits a 'current 
reversal' relative to the asymptotic direction of the circulation. The counter-circulating current in the core is due to 
a counter-moving bound state that appears below the Fermi level and dominates the current for distances of order 
< r 2£o- At high temperatures, T — > T c , this counter-moving bound state is thermally depopulated with the 
result that the current reversal in the core disappears in the Ginzburg-Landau limit. In summary, we find that the 
Andreev bound states dominate the current of vortices on the scale of a few coherence lengths. The nonequilibrium 
properties of vortices on this scale are expected to be dominated by the spectral evolution and dynamics of these 
bound states. 



ACKNOWLEDGEMENTS 

This work was initiated when the authors were participants at a workshop at the Institute for Scientific Interchange, 
Villa Gualino, Torino. The research of DW was supported in part by the Engineering and Physical Sciences Research 
Council, while that of DR and JAS was supported in part by NSF grant DMR 91-20000 through the Science and 
Technology Center for Superconductivity, the Max Planck Gesellschaft and the Alexander von Humboldt Stiftung. 
We also thank Dr. M. J. Graf for his comments on the manuscript. 



14 



IV. APPENDIX 



A. A matrix element 



In this appendix we derive a form for a matrix element used in section [ID. The matrix element in question is 
(C| [ii)fP£, S(e — Ha)]\C)i,i where we write 



We have 



F(f) 

H A = v f p c r 3 + A, A = A — ^(fiC + f 2 7 ? ). 



\iv f p c ,5(e - H A )} = i [f 3 (H A - A)S(e - H A ) - S(e - H A )(H A - A)f 3 



(44) 



(45) 



Then 



i[f 3 e, 5(e - H A )} - z{r 3 A, S(e - H A )} . 



(t\[iv fP( ,8(e-H A )}\Q 1 , 1 =iTr 



:(1 + r 3 )(C\([rse, S(e - H A )} - {f 3 A, 5(e - H A )})\Q 



(46) 



(47) 



Tr i({\r 3 A6(e-H A )\Q . 



Substituting the explicit form for A yields the relation 



(t}[iv f p ( ,5(e-H A )}\Q 1A = A Q ^Tr kr,- , m )\\M< - 7/ujO 



(48) 



(49) 



B. Approximation of bound states at large distances 
Equation ( |l2| ) gives the current density in terms of the Andreev Hamiltonian (pj) whose eigenvalue equation reads 

-ivfd^ + A - 



V( 2 + V 2 



V-(C) = Ei>(Q. 



(50) 



The parameter rj appearing in the above equation has the semiclassical interpretation as a c-number impact parameter. 
In this appendix we present approximations to the above equation for large values of the impact parameter, rj. 

For |?7| 3> £o we are justified to replace F(y/ ( 2 + rj 2 ) by its asymptotic value of unity. Furthermore, making the 
somewhat crude approximation 



m 9 ~> n-j-r + r 2 sign (rj) 



yields for the bound state 



-00 (C) — constant x exp I — t-tC 2 

• 2vf\n\ 



V2 / 



E ~ -A sign(77). 



(51) 



(52) 



This indicates that at large impact parameters, the bound states of the Andreev equation are found very close to the 
threshold of the continuum. 

When required, more refined estimates may be obtained by writing 



(riC + hy) 



T\ exp [i arctan(77/C)f 3 ] 



(53) 
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and then performing a unitary transformation Ha — * UHaU 1 , ip — ► tp = with ?7 = exp [| arctan(r;/£)f3] to 
remove the phase from the order parameter. The transformed Andreev equation is 



2 C 2 + V 2 



A fi 



^(0 = EtP(Q . 



(54) 



This is a one dimensional Dirac equation with a weak scalar potential, which has weakly bound states with energies 
near ±A . A "non-relativistic" treatment is appropriate in this case and we approximate the Dirac equation by a 
Schrodingcr equation. For example, for rj < we write 



ip = ipi 



1 



(55) 



with tpL,s scalars. Straightforward manipulations indicate that tpL approximately obeys the Schrodingcr equation 



v ) d 2 v f \y\ 
2A d( 2 2 C 2 + V 2 



MQ = {E-Ao)MO- 



(56) 



All the machinery of Schrodinger theory may be used on this equation to estimate e.g. the bound states. We can put 
a lower limit on the bound state energy. This may be obtained from the fact that the eigenvalues of the Schrodingcr 



operator are > V m i n , the minimum of the potential. Thus, E — A > Min^ 

v f 1 



v j \v\ 

2 C+n' 2 



i.e. 



E > A - 



2 \n\' 



T] < 0. 



(57) 
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